emg$ID = sub('.*(\\d{3}).*', '\\1', emg$Subject)
emg_sub =emg %>% select(Subject, ID, project, Stimulus,TAmplitudes,NumberOfTimesPresented, StimulusNumber) %>% unique() # to be added  Q_LTE_di, Q_CTQ_di, Q_STAIT_sum, Q_BDI_sum
emg_sub$picNum<-str_sub(emg_sub$Stimulus,11,-1L)
emg_sub$Stimulus = NULL
emg_sub$ID = gsub("[^0-9.-]", "",emg_sub$ID)
emg_sub$ID = gsub("^0", "", emg_sub$ID) 
emg_sub$ID = gsub("^0", "", emg_sub$ID) 

# prep df_rat
df_rat$picNum<-str_sub(df_rat$picture,1,6)
df_rat$picture = NULL
df_rat_emg = merge(emg_sub, df_rat, by = c("ID", "project", "picNum"))
#detach(package:plyr) # if "sanity" is empty
sanity = df_rat_emg %>% group_by(Subject) %>% tally()

dfEMG<-df_rat_emg%>%select(ID,project, picNum,TAmplitudes,NumberOfTimesPresented,Category, StimulusNumber,arousal_rating, valence_rating)

# get old ID structure back, if needed 
dfEMG$ID <- paste(dfEMG$project ,dfEMG$ID ,  sep = "_")
dfEMG$project = NULL 
emg_rat = dfEMG
# new: n = 401
hist(as.numeric(dfEMG$arousal_rating))

hist(as.numeric(dfEMG$valence_rating))

Outline

RSA from Ventura-Bort et al., 2022:, data from Koppold et al., 2023

Affective Startle Modulation Task

STR Based on T Amplitudes

Model Testing of Valence: Nearest Neighborhood, Anna Karenina (AK) and inverted AK

Model Testing Beyond Valence

Startle based on T Amplitudes

STR_Valence_RSM

Models of Valence

the RSMs of SCR will be tested against two models:

the Nearest Neighborhood model which assumes that the closer trials are in valence, the more similar they are
from the Anna Karenina (AK) model which assumes that trials with high valence will be more similar to each other than those rates as low arousal, we will derive the inverted model (Inverted AK Model), which assumes that lower values of valence are similar to each other, but higher values are more dissimilar

Behav_Single_STR_Valence_RSM

Behav_Single_STR_AK_Valence_RSM

Behav_Single_STR_invAK_Valence_RSM

Model tested againt each Model

ValRatingsNN_STRSingle_IndividualPlot

          ValRatingsAK_STRSingle_IndividualPlot

          ValRatingsinvAK_STRSingle_IndividualPlot

ggarrange( STR_ValRatingsNN_PermOutput[["Output"]][["PermutationPlot"]],
           STR_ValRatingsAK_PermOutput[["Output"]][["PermutationPlot"]],
           STR_ValRatingsinvAK_PermOutput[["Output"]][["PermutationPlot"]],
           labels = c("NN", "AK", "invAK" ),   hjust=-2,
           ncol = 1, nrow =3, widths = c(10,10))

Test the unique contribution of each model (i.e., NN and AK)

VECTOR

ggarrange(STR_Single_ValRatings_AK_ValRatings_NN_IndividualPlotRegression[[1]],
          STR_Single_ValRatings_AK_ValRatings_NN_IndividualPlotRegression[[2]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

ggarrange(STR_Single_ValRatings_AK_ValRatings_NN_PermOutputRegr[["Ouput"]][["RSMPlotRegressionM1"]],
          STR_Single_ValRatings_AK_ValRatings_NN_PermOutputRegr[["Ouput"]][["RSMPlotRegressionM2"]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

ggarrange( STR_Single_ValRatings_AK_ValRatings_NN_PermOutputRegr[["Ouput"]][["PermutationPlotM1"]],
           STR_Single_ValRatings_AK_ValRatings_NN_PermOutputRegr[["Ouput"]][["PermutationPlotM2"]],
           labels = c("iAK", "NN"),   vjust=(12),
           ncol = 1, nrow =2, widths = c(5,5))

Relationship between Similarity measures of the NN and invAK Model

STRCorrModComp

### Plot individual participants as example:

High in AK and and in NN Model

#Plot VP 
#-------------- Plot RSM Matrix
dataOriginal<-STR_corr_prep_all#data coming up from RSM functions it should contain the variable mean_all
mypalette<- numeric()#colors
minimum<-numeric()#scale
maximum<-numeric()# scale
Category<- "trials" # VP or trials
Participant<- "VP3Z02_1" # factor() #if individual participants are to be selected,
sortType<- "Valence" #variable to sort on
PlotMatrix(dataOriginal, minimum,maximum,Category,sortType,Participant, mypalette)
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(Participant)
## 
##   # Now:
##   data %>% select(all_of(Participant))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

High in invAK and close to 0 in NN Model

#Plot VP 
dataOriginal<-STR_corr_prep_all#data coming up from RSM functions it should contain the variable mean_all
mypalette<- numeric()#colors
minimum<-numeric()#scale
maximum<-numeric()# scale
Category<- "trials" # VP or trials
Participant<- "VP3Z02_68" # factor() #if individual participants are to be selected,
sortType<- "Valence" #variable to sort on
PlotMatrix(dataOriginal, minimum,maximum,Category,sortType,Participant, mypalette)

Low in invAK and moderate in NN Model

#Plot VP 
dataOriginal<-STR_corr_prep_all#data coming up from RSM functions it should contain the variable mean_all
mypalette<- numeric()#colors
minimum<-numeric()#scale
maximum<-numeric()# scale
Category<- "trials" # VP or trials
Participant<- "VP3Z02_18" # factor() #if individual participants are to be selected,
sortType<- "Valence" #variable to sort on
PlotMatrix(dataOriginal, minimum,maximum,Category,sortType,Participant, mypalette)

Model Testing Beyond Valence Is the AK the best fitting model?

Model Time

STR_TimeSortedByValence_IndividualPlot

Test whether the AK model is contributing beyond time

SCRSingleCorrModCompTime

ggarrange(STR_Val_AK_Time_NN_IndividualPlotRegression[[1]],
          STR_Val_AK_Time_NN_IndividualPlotRegression[[2]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

Model Category (Pleasant, Unpleasant, Neutral)

STRSingle_Categorysortedbyvalence_IndividualPlot

Test whether the AK model is contributing beyond Category

SCRSingleCorrModCompCat

ggarrange(STRSingle_Valence_AK_Category_Factor_IndividualPlotRegression[[1]],
          STRSingle_Valence_AK_Category_Factor_IndividualPlotRegression[[2]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

Model Valence

STRSingle_AroSortedByValence_IndividualPlot

Test whether the AK model is contributing beyond Valence

SCRSingleCorrModCompVal

ggarrange(STR_Valence_AK_Arousal_NN_IndividualPlotRegression[[1]],
          STR_Valence_AK_Arousal_NN_IndividualPlotRegression[[2]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

Model Category based on Arousal (Arousing, non-arousing images)

STR_Category2SortedbyValence_IndividualPlot

### Test whether the AK model is contributing beyond Category 2

SCRSingleCorrModCompCat2

ggarrange(STR_Valence_AK_Category2_Factor_IndividualPlotRegression[[1]],
          STR_Valence_AK_Category2_Factor_IndividualPlotRegression[[2]],
          labels = c("", "" ),   vjust=-7,
          ncol = 2, nrow =1, widths = c(5,5))

save.image("RSA_STR_Sound_VB22_Def.RData")
library(ggridges) # ridgeline plot
library(tidyr) # long format
library(colorBlindness)
library(RColorBrewer)
#install.packages("cols4all", dependencies = TRUE)
#library(cols4all)

long_data <- gather(SingleSimComparisonCat2, key = "model", value = "Correlation", 2:3)
ggplot(long_data, aes(x=Correlation, y=model))  + geom_density_ridges()
## Picking joint bandwidth of 0.0237

### single model alternative based on RSM data 
#STR_VariableModel1_VariableModel2_tableRegression = tableCorrelation
# long_data2 <- gather(STR_VariableModel1_VariableModel2_tableRegression, key = "model", value = "Similarity", 2:3)
# ggplot(long_data2, aes(x=Similarity, y=model))  + geom_density_ridges()+
#   scale_y_discrete(labels = c(VariableModel1, VariableModel2))


### all models
# test = cbind(STR_Val_AK_Time_NN_tableRegression,
 # STR_Valence_AK_Arousal_NN_tableRegression,
# STR_Valence_AK_Category2_Factor_tableRegression, 
# STR_Single_ValRatings_AK_ValRatings_NN_tableRegression,
# STRSingle_Valence_AK_Category_Factor_tableRegression)

# test2 = test[-c(4,7,10,13)]

test = cbind(STR_Val_AK_Time_NN_tableRegression,
        STRSingle_Valence_AK_Category_Factor_tableRegression,
        STR_Valence_AK_Category2_Factor_tableRegression, 
        STR_Valence_AK_Arousal_NN_tableRegression)
#test2 = test[c(1,3,6,9,12)]#
test2 = test[c(1,2,5,8, 11 )]
names(test2)<-c("VP", "Time", "Cat", "Cat2", "Arousal")

long_data2 <- gather(test2, key = "regressor", value = "rho value", 2:5)
ggplot(long_data2, aes(x= `rho value`, y=regressor, fill = regressor))  + geom_density_ridges()+
  geom_density_ridges(quantile_lines=TRUE,
                      quantile_fun=function(x,...)mean(x), color = "white") +
    scale_fill_brewer(palette = 4) +
  theme_ridges() + theme(legend.position = "none")+
  geom_vline(xintercept = 0, color = "red",size = 1.5, linetype = "dashed")
## Picking joint bandwidth of 0.0296
## Picking joint bandwidth of 0.0296
## Warning: Using the `size` aesthietic with geom_segment was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# colorblind friendly
cols = c("#E69F00"  ,   "#004949"  ,   "#009292" ,    "#FF6DB6"    , "#FFB6DB" ,    "#490092"  ,   "#006DDB" ,
    "#B66DFF"  ,   "#6DB6FF"   ,  "#B6DBFF")
ggplot(long_data2, aes(x= `rho value`, y=regressor, fill = regressor))  + geom_density_ridges()+  geom_density_ridges(quantile_lines=TRUE,
                      quantile_fun=function(x,...)mean(x), color = "black") +
  theme_ridges() + theme(legend.position = "none")+
    scale_fill_manual(values = cols)+
  geom_vline(xintercept = 0, color = "red", size = 1.5, linetype = "dashed")
## Picking joint bandwidth of 0.0296
## Picking joint bandwidth of 0.0296

#displayAvailablePalette(color="white")
library(dendextend)
## 
## ---------------------
## Welcome to dendextend version 1.16.0
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------
## 
## Attaching package: 'dendextend'
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## The following object is masked from 'package:stats':
## 
##     cutree
library(circlize)
## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================
my_data = Behav_Single_STR_AK_corr_prep_all
# Select the columns for clustering
data_for_clustering <- my_data[, c("var1", "var2", "mean_all")]

# Perform hierarchical clustering
dist_matrix <- dist(data_for_clustering)  # Calculate the distance matrix
hc <- hclust(dist_matrix)  # Perform hierarchical clustering

# Plot the dendrogram
#plot(hc, main = "Dendrogram of Clustering")

# Convert hierarchical clustering object to a dendrogram
dendro <- as.dendrogram(hc)

dendro <- dendro %>% 
  color_branches(k=4) %>% 
  color_labels

plot(dendro)

# Create a circular dendrogram
# plot the radial plot
par(mar = rep(0,4))
# circlize_dendrogram(dend, dend_track_height = 0.8) 
circlize_dendrogram(dendro, labels_track_height = NA, dend_track_height = .4)